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1. INTRODUCTION 


The purpose of this work is to investigate stability of the spinning 
Skylab control system [1, 2] by the modem decomposition-aggregation methods 
[3, 4]. Such an investigation is motivated by the fact that the mathemati- 
cal model of the system is of high dimension and a straightforward analysis 
would become bogged down in the welter of detail requiring an excessive com- 
puter storage and time to complete the investigation. The multi-level decom- 
position-aggregation approach offers to solve the stability problems "piece- , 
by-piece" and not only make more economical the computer use, but also reduce 
the liability of the errors in the analysis. Furthermore, by decomposing the 
system into parts that have important physical meaning, the decomposition- 
aggregation approach yields significant structural information about the be- 
havior of the system, which is not generally available in a straightforward 
stability investigation. 

The outline of the work is divided into two parts: Theory and Applica- 
tion. The part on th&ofiij presents the mathematical basis of the decomposition- 
aggregation approach - the concept of vector Liapunov functions [3, 4]. In 
Section 2, it is shown how a vector differential equation of high dimension,, 
which describes a large dynamic system, can be decomposed into a number of 
vector differential equations of lower dimensions, which represent intercon- 
nected subsystems. Then, in Section 3, a method is outlined by which the 
stability properties of each subsystem are aggregated into a single Liapunov 
function. The vector Liapunov function is formed which has subsystem Liapunov 
functions as components. A linear vector differential inequality is then 
formed in terms of the vector Liapunov function, which represents the aggre- 
gate system model. The matrix of the model, which reflects both the stability 
properties of the subsystems and the nature of their interconnection is ana- 
lyzed in Section 4 to conclude stability of the over-all system. 

i 
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The appLLcation part of the outline presents an application of the 
decomposition-aggregation method to determine stability of the spinning 
Sky lab control system. In Section 5, equations of motion [1, 2] of the 
system are given which include both the passive stabilization by extend- 
able booms with tip masses and the active stabilization by control torques 
1 about the body fixed axes. Stability analysis of the passive control in 
' Section 6, starts by decomposing the equations of motion into two sets of 
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equations which describe the wobble motion and the in-plane motion. Then 
two sets of equations are treated as subsystems which makes the coordinates 
of the tip masses to appear explicitly in the interconnections as structur- 
ally important coupling parameters. The decomposition-aggregation method 
is then used to determine stability of the over -all system as a function of 
the coupling parameters. The entire stability procedure can be conveniently 
programmed on the digital computer. The programs and their description is 
given in the Appendix. In Section 7 the decomposition- aggregation method 
is used to determine stability of the Skylab control systems when control 
torques are used. The torques are considered as linear functions of the 
.states and are applied about the corresponding body fixed axes. The unspeci- 
fied parameters of the linear functions provide a considerable freedom which 
can be used in the stabilization of the vehicle on the subsystem level. The 
task of establishing stability of the entire vehicle is, therefore, increas- 
ingly easier than that in passive control. However, how to intelligently use 
the available freedom in the control of the subsystems and produce a higher 
stability degree of the over-all system, is not yet satisfactorily resolved. 

A solution to this problem represents one of the major goals of the future 
efforts in development and application of the decomposition-aggregation 
methods for stability analysis of complex attitude control systems. 
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PART I 


THEORY 



2. DECOMPOSITION 


Decomposition of dynamic systems can be used to overcome the "large 
size" of the stability problems. A "large" problem is considered as con- 
stituted of coupled "smaller" problems which can be (when isolated) solved 
efficiently in sequence. Then, the solutions of the "smaller" problems 
are combined together with the constraints on couplings into an aggregated 
model which is relatively simple to solve. Decomposition algorithms pro- 
duce considerable saving in both computer storage and the time required to 
complete the solution of the original problem. Furthermore, if decomposi- 
tion is performed so that the "small" problems can be physically interpreted 
(e.g. motion in the pitch plane), a decomposition analysis may produce impor- 
tant structural information about the "large" over-all problem [1, 2]. 

Let us consider a continuous dynamic system S described by the vector 
* 

differential equation 

x = f(t, x) , (2.1) 

where x(t) £ if 1 is the state of the system S . The function f : T x R 
-*■ R 11 satisfies a global Lipschitz condition so that the solutions x(t; tg, 
Xq) of (2.1) exist and are unique and continuous for all initial conditions 
(tg, Xq) G T x if 1 and t G Tg . The symbol T stands for the time inter- 
val (x, + <*>) , where x is a number or the symbol - » , and 7"g is the 
semi-infinite time interval [t Q , + °°) . 


* 

With some obvious exceptions , lover case Roman letters denote vectors , capi- 
tal Roman letters denote matrices, and Greek letters denote scalars. 
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We assume that the function f(t, x) in (2.1) satisfies the condition 

f(t, 0) = 0 , Vt 6 T (2.2) 

and that the origin x = 0 of the state space R 11 is the unique equilibrium 
state of the system S . In the following development, the emphasis will be 
on global stability of the equilibrium x = 0 and the uniqueness of the equi- 
librium causes no reduction of generality in the results. 

A crucial assumption in the following stability analysis is that the dy- 
namic system S can be (conveniently) decomposed into s dynamic subsystems 
and described by the vector differential equations 

*i = x i^ + h i^ t} *) > i = 2, ... , s . (2.3) 

The free (uncoupled) dynamic systems are described by 

x = g i (t, x.) , i = 1, 2, ... , s (2.4) 

n. 

where x^(t) 6 R 1 is the state of . Therefore (2.3) describes the mo- 
tion of the interconnected (forced) subsystems S. , where the function g. : 
n. n. 1 1 

T x R •+■ R corresponds to the subsystem S. itself, and the function 
n. 1 

lu : T x T -> R 1 represents the action of the composite system S on the 
subsystem . We assume again that 

g i (t, 0) = 0 , Vt G T , Vi = 1, 2, ... , s (2.5) 

and = 0 is the unique equilibrium state of the subsystem 5^ . 

From (2.3) and (2.4), it is clear that the state x^(t) £ R 11 * of the 
subsystem is the i-th vzcton component of the state x(t) of the over- 
all system S . That is, x(t) can be written as 
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x(t) = [x[(t) X^(t) ... Xg(t)] T , (2.6) 

and the state space R n of the system S can be represented by the Cartesian 
product 

n. • n 9 n 

R 11 = R 1 x R z x ... x R s . (2.7) 

Furthermore, each state x^(t) of the subsystem S can be written as 

x i (t) = [x n (t) x i2 (t) ... x^ (t)] , * (2.8) 

where x.p ^2’ "' x in. are ' ica ^ z/l components of the vector- x^(t) . : 

To illustrate the decomposition formulated above, let us consider the 
motion of a disk fixed to a rotating shaft as shown on Fig. 2.1. The system 
is regarded as a massless elastic shaft with a mass particle attached at 
the center. Friction is assumed to be internal to the shaft. If p and y 
represent deflections of the mass particle in a coordinate system rotating 

at the angular velocity w of the shaft then the linearized equations of 

* , . . 

motion are * : - - - * ■ - 

. 2 ' • . : 

mp + bp + (c - mw )p - 2nui)y = 0 ... 

•• • 2 

my + by .+. (c„ - nui). ) y + 2mwp •?. 0 , . - (2.9).- 

where m is the mass, b is the damping coefficient, and c is the stiff- 
ness coefficient of the shaft. .• .... . 

Equations (2.9) describing the system S of Fig. 2.1, can be given in 
the state formas ' ' ' 

Ziegler, H. , "Principles of Structural Stability", Blaisdell, Waltham, Mass., 

1968. .. ... . ••• .■ .r • ■ . 
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In (2,11), 


* 1 (t) 


[x n (t) x 12 


(t)] T , 


x 2 (t) = [x 21 (t) x 22 (t)) 


( 2 . 12 ) 


are the states of the two subsystems and S 2 , which are (vector) compo- 
nents of the state vector 


x(t) = [x n (t) x 12 (t) x 21 (t) x 22 (t)] T = [xj(t) x 2 (t)] (2.13) 


for the entire system S . 

By comparing equations (2.11) with equations (2.3), we see that the free 
subsystems are described by 



(2.14) 


Interactions among the subsystems are represented by the functions 



0 

0 


0 

0 

^(t, x) = 

0 

2oj 

x 2 h 2 (t, x) = - 

0 

2ui 


(2.15) 
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Since the system S is decomposed into two subsystems and S 2 with 
apparent physical interpretation, a subsequent stability analysis can yield 
structural information about the system. 
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3. AGGREGATION 


The Liapunov direct method can be viewed as a process by which the stabi- 
lity properties of the motion in the vector space are aggregated into a single 
scalar function - the Liapunov function. The aggregation procedure, therefore, 
produces an essential reduction in dimensionality of stability problems at the 
price of a reduced accuracy in the results. It is also true that generation 
of appropriate Liapunov functions can seldom be automatic and becomes increas- 
ingly difficult for systems of high -dimens ion. Consequently, in stability ana- 
lysis of large-scale systems, a straightforward attack by a single Liapunov 
function becomes cumbersome and another level of aggregation is desirable. 

After a large-scale system is decomposed into a number of interconnected 
subsystems, a Liapunov function is assigned to each isolated subsystem. The 
subsystem Liapunov functions are then used as components of a vector Liapunov 
function to construct the aggregate model of the over-all system. The aggre- 
gate model is a linear (vector) differential inequality written in terms of the 
vector Liapunov function, which can be effectively examined for stability by 
known methods. 

Since the subsystems are of relatively low order and often with similar 
characteristics, the generation of appropriate subsystem Liapunov functions is 
not exceedingly complex. Moreover, the decomposition of the system may be di- 
rected towards producing the subsystems for which Liapunov functions with de- 
sired properties are available. 

It is, however, true that the necessary properties of subsystem Liapunov 
functions for construction of the aggregate model, require stability of the sub- 
systems Which constrains the decomposition-aggregation stability analysis. 
Furthermore, in forming the aggregate model, approximations are involved which 
increase the conservativeness of the over-all results beyong the extent intro- 
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duced on the subsystem level. 

StnblLiXy pAopeAtteA ofi the AubAyAtemi are aggregated by scalar func- 
tions v.: T x I! 1 + such that v^(t, x^) £ ' (T x R 1 ) , which have 

the following "linear" estimates: 


ri i 1 ll x ill - v i — n 12 1 l x i I 
; ’i ? -"usihH 

I I grad v ± \ I < n i4 , (3.1) 

T l/° 

where n. » ’ n i3 ’ n i4 are P os; '- t i ve numbers and | |x^| | = (x^x_^) 

is the Euclidean norm of the vector x^ . In (3.1 ) ; v . ' is the total time 
derivative of the function v^(t, x^) along the motion of the free subsystem 
described by (2.4), that is 


v i = 


It V i + (grad V i )T 


(3.2) 


Liapunov functions v^(t, x^) with estimates (3.1) are available for 
large classes of dynamic systems. For example, consider the class of systems 
with linear, time -invariant subsystems, with 

g i (t, x i ) = P i X;L , (3.3) 


where the P.'s are constant-coefficient, n. * n. matrices. Then, for each 

l ’ll ’ 

subsystem there exists a Liapunov function of the form [3] 


v H i x i > 1/2 


(3.4) 


whose time derivative on the trajectories of the uncoupled subsystem is given 
by 
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(3.5) 


v. = - (1/2) V : 1/2 x[ G. x. 


where " -G ± - th p^_ + P^^ H- , G-. > 0 


Direct calculations give 


(3.6) 


n il 

n i2 

n i3 

n i4 



1/2 

1/2 


(H t ) 

( H i) 


XfG^A' 172 ^.) 
1/2 (H i )A(H.) , 


(3.7) 


where A and A denote the minimum and maximum eigenvalues of the indicated 
matrices, respectively. ! 

In the following Section 4, we will establish the fact that the existence 
of a Liapunov function with estimates (3.1) implies and is implied by the ex- 
ponential stability property of the corresponding dynamic system. 

InteJULCJtlom, among tkz 6u.bAyAtm6 are assumed to satisfy the following 
"connical" constraints 

IlhiCt, x)||< I e.llx.11 (3.8) 

j=l J J 


where are nonnegative numbers. For example, if the interactions are 

linear, time invariant and of the form 

s 


h. (t, x) = y Q. .x. , 
1 j£l J 


(3.9) 


where the matrices (y are n. x n^ constant-coefficient matrices, then the 

£n £ 

1/2 


numbers can be taken as 


£• .= [A (Or .0. .)]■ 


(3.10) 
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To form the aggregate model of the system S using Liapunov functions 
v^(t, x^) with estimates (3.1) and interaction constraints (3.8), we take 
the total time derivative of v^ (t, x.) along the motion of the intercon- 
nected subsystem described by (2.3), that is, 

= v^ + (grad v.) T h^t, x) , i = 1, 2, ... , s . (3.11) 


By applying inequalities (3.1) and (3.8), we can rewrite (3.11) as 

■ -I ® -1 

v. < - n-o n -7 v. + n-/i y £•• v. , 

~i — i2 i3 i i4 jl j ’ 

i = 1, 2, ... , s . (3.12) 

Now, we can define an s vector v function v: T x R 11 R^ using 
(3.5) , 

v = (y 1 v 2 ... v s ) T (3.13) 

which has as components the Liapunov functions v^ related to each subsystem 
. The vector function v(t, x) , is called the vecloh. Liapunov function 
[l , 2]* With the notation (3.13) , scalar inequalities (3 .12) can be rewritten 
in a vector form as 


v < Av 


(3.14) 


where the s x s matrix A = (a^) has the elements specified by 


o - * - 1 * r -1 

a. • • “ ** 0 • • T| • T) • <7 * c • « H • -i T1 • a • 

ij ij i2 i3 s i;j l jl ‘i4 * 


(3.15) 


where 5^ is the Kronecker delta. 

Inequality (3.14) is a vector differential inequality and is referred to as 
the aggregate model of the system S . The matrix A is called the aggregate 
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maXxix. 

The aggregate model contains the necessary information about the stabil- 
ity properties of the system S . The dimension of the model is s which is 
less, or at most equal to the dimension n of the original system. This pro- 
duces the desired reduction of dimensionality in the stability problems assoc- 
iated with large-scale systems. 

For the example, taking jjj- = = 1 and oi = 0.04 we have 


P. 

1 


0 


- 1.0 



1.0 


i - 1, 2, 



0 0 
0 0.08 



0 0 
0 -0.08 


Then, the choice 



gives, through the application of (3.7) 


(3.16) 


' (3.17) 


- n il = 
n i2 = 
n i3 = 
n i4 = 


1.17 

1.90 

0.26 

1.63 


i = 1, 2 


(3.18) 


From (3.10) the numbers 



are found to be 
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(3.19) 


5 12 = hi 


0.08 . 


Using (3.15), the elements of A are calculated as 


a ll = a 22 ’ -°- 14 


a 12 “ a 21 " 0,11 ’ 


and the aggregate model (3.14) becomes 


(3.20) 


-0.14 

0.11 

0.11 

-0.14 


(3.21) 


How the aggregate model (3.21) can be used to conclude (exponential) sta- 
bility of the original system (2.10), will be explained in the following section. 
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4. STABILITY 


The purpose of this section is to show how stability of a large-scale 
system can be determined from the stability of the subsystems and the nature 
of their interactions. On the subsystem level, we first conclude that 
Liapunov function estimates (3.1) imply and are implied by exponential stabi- 
lity of the subsystems. Then, we establish the same kind of stability for 
the over-all system by demonstrating stability of the comparison (linear) sys- 
tem. 

Let us consider a free dynamic system described by the differential equa- 
tion 

*i = g i^’ ’ C 4 * 1 ) 

which represents one of the subsystems of Section 2. 

We use the following standard definition of exponential stability [ ] : 

definition 4.1 . The equilibrium state x^ = 0 of the free subsystem Is 
globally exponentially stable. If and only If there exist Two positive, numbers 
ou and 3^ Independent o f the Initial conditions (tg, Xg) such that 

I I x ± (t: ; t Q , Xg) 1 1 < ou||x i0 || exp[-3.(t - t Q )] , 

n. 

Vt G Tg , V(tg, x i0 ) G T X R 1 (4.2) 

Exponential stability can be established by the following modification 
[ 4 ] of the well-known Krassovskij [ 5 ] result: 

Theorem 4.1 . The equilibrium stale x^ = 0 of the knee subsystem Is 

globally exponentially stable If and only If there exists a function v. (t, x.) 

n. 1 1 

on T x R 1 with estimates (3.1). 

This Theorem follows directly from Krassovskij ' s result when one con- 
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siders vV 2 . instead of v^ . This modification, however, allows a simple 
construction of the aggregate model as shown in the preceeding Section 3 . 

Let us prove the sufficiency part (the "if" part) of Theorem 4.1. From 
the estimates (3.1), it is easy to write the following scalar differential in- 
equality 


v. 

l 



i3 


v. 

l 


(4.3) 


Inequality (4.3) can be integrated to yield 


V*. V t; V x 10 }1 -W x i0 ] ex Pl-"i2 

n. 

Vt 6 T q , V(t Q , x. Q ) 6 T x R 1 . (4.4) 


By applying again estimates (3.1), we can rewrite (4.4) in terms of 
as 

.■1 _ II.. II r -1 


1^1 


l x .(t; t Q , x i || < 


n il n i2 H x iO^ ex P[' n i2 n i3 ( - t ' t 0^ 


n. 


Vt G T q , V(t Q , x iQ ) 6 T x R 1 . 


(4.5) 


Comparing (4.5) with (4.2), we get the positive numbers a and & as 


a i = n ii n i2 


B i = n i2 n i3 


(4.6) 


which are independent of initial conditions (tg, x^g) . This proves the suf- 
ficiency part of Theorem 4.1. 

Let us now consider the system S described by (2.1) and its aggregate 
model given by (3.7). We prove the following: 

Tk 2 . 0 A.rn 4.2 . The zquAZibsuum &tat& x = 0 the. 6y&tm S ' -Lt> globaZZy ex- 
ponenttaJtly stable, the. s * s aggregate mouUux A = (a^) AatiA&ieA thz 
tne.quaLitieA - 
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(4.7) 




a ll 

a l2 

a ls 

a ll a l2 

> o, ... , (-1) S 

a 21 

a 22 

**• a 2s 

a 21 a 22 


a sl 

a s2 

a 

ss 


psioot . According to Definition 4.1, we should prove that inequalities (4.7) 
represent sufficient conditions for the existence of two positive numbers a 
and g independent of the initial conditions (tg, x Q ) , such that the solu- 
tion x(t; tg, Xq) of equation (2.1) satisfies the following inequality; 

i l x (t; t Q , x 0 )|| la||x 0 || exp[-g(t-tg)] , Vt G T Q (4.8) 

for all (t Q , Xg) 6 T * R 11 . 

Let us consider a decrescent, positive definite, and radially unbounded 
function v: R s , 

v(v) = b T v , (4.9) 

where b > 0 is a constant s -vector, as a candidate for Liapunov's function 
[9] for the system S . 

Taking the total time derivative along the solutions x(t; tg, x Q ) of 
(2.1), we get 

v(v) = b T v[x(t; t Q , Xg)] .. . . (4.10) 

. . T * ... 

Premultiplying inequality (3.7) by b > 0 , we obtain 

C ^b T Av . . . •- (4.11) 

Let us rewrite (4.11) as 
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s s 

I \ l 
i-1 j=l 

s s 


v < j b. y a. . v. 
— j“i i ij J 


< y v. y b. a.. 
-j=l J i=l 1 13 


1 -i v j b jl a iii + l v j l b i a.. 
j=l 3 3 33 j=l 3 i=l 1 13 

i£j 


(4.12) 


From the definition (3 .15) of coefficients 
elude that 



in the matrix A 


we con 



< 0 i = j 


1 0 i t j , 

(4.13) 


that is, the matrix A has nonnegative off-diagonal elements. This special 
structure of A makes the conditions (4.7) necessary and sufficient [6] for 
the existence of a vector b > 0 > 0 , Vi = 1, 2, ... , s) and a number 

$ > 0 such that 

-I s 

|a-. | - b. I b. a.. >_ g , Vj = 1, 2, , s (4.14) 

33 3 i=l. 1 13 

From (4.12) and (4.14), we get 

v < Bv , Vt 6 T q , Vv 6 R* . (4.15) 

Integrating inequality (4.15), we obtain 

v[x(t; t Q , x Q )] <v(x 0 ) exp[-3(t-t 0 )] , Vt G T Q , Vv 6 rJ . (4.16) 

It is left to show that inequality (4.16) implies inequality (4.8) and, 
thus, global exponential stability of S . The left-hand side of inequality 
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(4.16) can be developed as 


v[x(t; t Q , x Q )] = ^ v.. 


— b m "ml X IM 
1=1 

i b m "ml IM V X 0 J H • 


(4.17) 


where 


b m = T b i > \a = “ n "u . 


(4.18) 


In developing (4.17), use is made of the estimates (3.1) and the fact that 

I IM £ f 

1=1 1 i=l 1 

written on the basis of the following derivation 


2 2 

x - 1 1 ~ ! |x| | . The right-hand side of (4.16) can be re- 


v[x 0 ] - b t V. 


= l b, |v, 


i=l 


l ' l ' 


lbw Ji ,V i ! 

i sl/2 ^jl|v|| 

- S ^ ^ ^2 1 ’ 


(4.19) 


where 


= max ^ , ^2 = max n i2 


(4.20) 


In derivation (4.19), again the estimates (3.1) are used and the well-known 
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relationship concerning the Euclidean norm and the absolute value norm, 
s 1/2 | |v| | > |v| . 

By using (4.16) and (4.18), we can rewrite inequality (4.15) as inequa- 
lity (4.8) with 


1/2 . , -1 -1 
° ■ 5 Vm n M2 n ml 


(4.21) 


and 6 defined in (4.13). This proves Theorem 4.2. 

To calculate a and 8 in (4.8), we need first to determine a -vector 
b > 0 in (4.9) , 

v = b*v , (4.9) 

and. (4.11) which we can rewrite as 

v <_ - c^v (4.22) 

where the vector c is defined as 

c T = -b T A . (4.23) 

Since the coefficients a^ satisfy (4.11), conditions (4.7) are neces- 
sary and sufficient [ 6 ] for existence of a positive vector b for any posi- 
tive vector b for any positive vector c . Consequently, one can calculate 

a vector b from 


b 


T 


T 

-cA 



(4.24) 


Now, the .constant a is calculated from (4.21) and the constant 8 is 
calculated from 


8 = min { | a*. 


-b: 1 

3 


I b. a. . } 
i=l 1 ^ 


(4.25) 
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which is equivalent to (4.14). 

Let us consider again the specific example of Section 3. First, to con- 
clude global exponential stability of the system, we verify that the aggregate 
matrix 

“0.14 0.11“ 

A = (4.26) 

_0.11 -0.14_ 

of (3.21) satisfies conditions (4.7). . 

By choosing r . .. 

c T = [1 1] , (4.27) 

and using (4.24), we calculate 

b T = -33.4 [11]. (4.28) 

With this b , we compute a and 8 from (4.21) and (4.25) as 

a = 2.30 , 0 = 0.03 . (4.29)’ 


23 



PART II 
APPLICATIONS 



5. SKYLAB MODEL 


The Sky lab [ 1 ] as an earth-orbiting manned space station, which is 
designed for prolonged space flights, may be required to provide artificial 
gravity environment. Therefore, Marshall Space Flight Center of NASA initi- 
ated the feasibility of spinning the Skylab about a principal axis of inter- 
mediate moment of inertia and producing the artificial gravity effect [ 1 ] . 
Since such spin cannot be achieved without stabilization, it was hoped that 
passive stability could be established by deploying masses either on cables 
or extendable booms as shown on Fig. 5.1. Such configuration has the princi- 
pal axis of maximum moment of inertia pointing (in the same direction as the 
solar panels) to the sun. 



Fig. 5.1 
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A simplified model [ 1 ] of the spinning Skylab consists of a core mass 
with two tip masses connected to it by flexible massless beams lying in two 
different planes as shown on Fig. 5.2. 



The angular velocity vector of the vehicle may be written in body-fixed 
coordinates 1, 2, 3, as 

t w l w 2 w 3 + ft] T , (5.1) 

where |w^| « 1 (i = 1, 2, 3) represent small perturbations about the steady 

state velocity ft . Small displacements of the two tip masses m from the 

k 

steady state are denoted by u.(i s 1, 2, 3; k - 1, 2) . The rotational dyna- 
mics of the Skylab may be represented by a set of nine differential equations 

if 

written in terms of w^ and u^ . It is possible to reduce the set of nine 
equations to six equations by using the substitution 



where u^ now represents the skew symmetric mode of the elastic deformation 
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and hence causes angular motion about the vehicle’s steady-state attitude. 
Since stability of rotational motion will be of interest, only the skew sym- 
metric mode is considered. 

The linearized equations of motion are 


f Vl + ^3 ‘ 1 2 )0m 2 + mr 2 ((i 3 + 0 U 3 } 


wobble < 
motion 


? 

-mr^ (2flu^ + U 2 ) = T^ 

(1^ - I 3 ) Ow^ + I 2 ^2 + mr 3 ^1 ~ U 1 " 2^2^ = T 2 


2mr 


2 (w^ + Ow 2 ) + ^3 + dj lij + (kj + m£r) u^ =0 


(5.3) 


r ^ • 


spin J 
motion ' 


X 3 W 3 " mr 2 ^ U 1 " 2 fiu 2^ T 3 
2mr 3 (Ow^ + v^) - 2mr 2 W2 + noi^ 


+ d i u i + k i u i " 2m«U2 “ 0 


2mr 3 (-w^ + «w 2 ) - 4mr 2 nw 3 + 2mRu- 1 + mu 2 


+ d2 U 2 + (T«2 " ) u 2 = 0 , 


(5.4) 


where: 1 ^, 1 ^, 1 ^ are the principal moments of inertia in the steady state 

( 1 ^ < I 2 < Ij) ; k^ are the stiffness coefficients of the nonrotating booms; 
mft is the geometric stiffness introduced by spin (the overall boom stiffness 
in the. 1, 2, 3 directions is, therefore, k^ + mP , k 2 , k^ + mft , respective- 
ly); d^ (i = . 1 , 2, .3) are damping coefficients relating the structural damp- 
ing to elastic deformation velocities; r. , -I\ (i = 1, 2, 3) are the coor- 
dinates of the two tip masses at the equilibrium, respectively; and (i = 

1, 2, 3) are the applied torques about the body -fixed axes. 
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Stability analysis of the Skylab model (5.3) will be performed separately 
for two controls, namely, paA&ive. control when all the torques are identically 
zero T. = 0 (i = 1, 2, 3) , and aejbLvz aont/iot, when the torques are linear 
functions of the states of the Skylab. 

In stability analysis, the following Nomenclature and Physical Character- 
istics of the Skylab are used: 

Nomenclature 

principal moment of inertia of body about i coor- 
2 ^ a 2 

dinate 1^ + Zmr^ , , Ij + 2 mT 2 , respectively 

principal moment of inertia of rigid core body about 
i body - fixed coordinate 

stiffness coefficient characterizing nonrotating 
boom stiffness 

m = tip mass of boom 

T\ = applied torque about i coordinate 

t = time 




( ) = d/dt = differentiation with respect to real time t 



m 

l 


= skew symmetric mode of elastic deformations 

= displacement of m tip mass from spinning steady 
state in i direction (m = 1, 2) 



= perturbation (about spinning steady state) velocity 
about i coordinate 
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r 


2 


r 


3 



1+K 1 T 2 


Y = 




u . = u./2r„ 

i i 2 


a? = k./mfi^ 
l l 



T = fit 


fi 


V. 

1 


w. 

1 


( ) = d/dx 


= steady-state boom dimension in 2-axis direction 
from center of mass to tip mass 

= the asymmetry in the setting of the’ booms 

= ratios of inertia (I^IjVl^ and ^3 -I l^ I 2 
respectively 

J 2 

= ratio of inertia y^- 

L 1 

= dimensionless inertia ratio 

= dimensionless damping ratio 

= general skew symmetric coordinate 

= dimensionless natural frequency coefficient of 
boom 

= dimensionless length ratio 

= dimensionless time 
= steady-state spin rate about 3 axis 
= dimensionless wobble ratio (i = 1, 2, 3) 

= differentiation with respect to x 


subscript i 


= index referring to three body-fixed coordinates 
(i = 1, 2, 3) 
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2 

Yj = Zmr^/Ij = dimensionless inertia ratio 
*3 

8 - f- = dimensionless inertia -ratio 
L 1 

Physical Characteristics 

^ - 1.25 x 10 6 kg m 2 
I 2 = 6.90 x 10 6 kg m 2 
I_ = 7.10 x 10 6 kg m 2 

T 1 = 0 
= 23.3m 

r 3 = -1.53m 

m = 227 kg 

k^ = k 3 = 146 N/m 

k 2 = 7.4 x 10^ N/m 

d 1 = d 3 = 0.04 (k 3 m) 1</2 

d 3 = 0.04 (k 2 m) 1/2 

n = o.6s"^ 


32 



6. PASSIVE CONTROL 


By adding the extendable booms with tip masses to the Skylab (Fig. 5.1), 
the spinning vehicle meets the condition 

h < h ■" h > 

and can be stably spun about the 3-axis. It is of interest in this section, 
to study the stability properties of the system model . 3) when the active 
control is not present, that is, 




An important feature of equations (5.3-4) and (6.4-5) is that when = 0 , 
they become uncoupled into two sets of equations : the wobble motion (w^ , w^ , u^) 
described by (5.3) or (6.4) ; and spin motion (w^, u^, u 2 ) described by (5.4) 
or (6.5). This suggests that the influence of the asymmetry in the arrangements 
of the booms (r^ f 0) can be treated as the coupling parameter between the two 
motions. In the decomposition-aggregation analysis, each motion represents a 
subsystem. 

Passive control equations (6.4-5) can be rewritten as follows: 

I I? f It 

V 1 " K 1 v 2 + y ^ p 3 + y 3^ ' + y 2 " u 2^ = 0 

I ft I 

- ^2 aV l + aV 2 + ^ Y ^ p l ’ y l ’ 2y 2^ = 0 

i n i 2 

+ v 2 + u 3 + A 3 u 3 + (ct 3 + 1) P 3 = 0 (6.6) 

i ii '2 * 

5(vi + \> 2 ) + p^ + A-jPj + - 2y 2 = 0 

^(-v^ + v 2 ) + 2p 1 + p 2 + A 2 2y 2 + (°2 ’ ~ 0 * (6.7) 

where the notation is introduced as in the Nomenclature. The dimensionless 
parameter 5 = r 3^ r 2 t ^ e C0U P^^ n S parameter between the two sets of equa- 
tions (6.6) and (6.7). 

The state space representation of the over-all system S described by 
(6.6-7), is obtained as 

S: x’ ( t) = P x (x) , (6.8) 

where the state 8-vector x(t) is chosen as 

! I t 'J' 

x (t) = (v-j^ v 2 P 3 P 3 v 1 P x P 2 P 2 ) , (6.9) 
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and the 8 * 8 matrix P is given in Fig. 6.1. 

In order to obtain the state equations (6.8), it is necessary to write 
equations (6.6-7) in the vector form 

Bx (t) + Cx(x) = 0 (6.10) 

where x(x) is the state vector (6.9), and the 8x8 matrices B and C 
are given by 




The matrix P of (6.8) and Fig. 6.1 is obtained from (6.10) as 

P - - B _1 C . (6.12) 

The system S of equation (6.8) can be decomposed into two interconnected 
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Fig. 6.1 



subsystems described by 

Sta S l : x i« ■ VlW * £ V» x i« ♦ 5Q 12 (0 * 2 (t) (6.13) 

motion V x 2« = W t) * N 2 1 (£) x l« * £ V« x 2« (6 - 14) 

where the state vectors x(t), x^(t), X 2 (t) of the system S and the two sub- 
systems and S 2 are 

“*!«] X 1 (T ^ = v 2 y 3 y 3 )T 

X(t) = ; , I rp 

" x 2^ t -* x 2^ t ^ = ^ y l y l u 2 y 2^ * (6.15) 


In (6.13-14) j the 4 x 4 matrices and P 2 correspond to the subsystems 

S-p and S 2 , and the 4x4 matrices Q n (0, Q 12 (0, Q 22 te)* Q 21 (0 repre- 
sent the interconnections between the two subsystems : 



n t, 1 r, 1 J- 

0 p 12 p 13 p 14 

P2! 0 00 

0 0 0 1 

n 1 1 1 

0 p 42 P 43 . P 44 


P 12 = ' P 42 = (K 1 + Y)/(1 * y) 
P 13 = y °3 /(1 ' y) i 
p 14 = A 3 y/(1 " y) 

P 21 = K 2 

P 43 = ‘ f g 3 + 1 ‘ Y)/(l- y) 

P44 = ‘ A 3 /(1 ‘ Y) 


"0 10 0 “1 p 21 = - a 2 

2 2 2 

P 21 P 22 0 P 24 P 22 = - A 1 

0... 0 0 1 P | 4 = 2 

_ 0 P 42 P 43 P 44_ p 42 = " 2 

p 43 = -(*2 ‘ 15 
P 44 A 2 
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11 _ 11 ( 2y + K 1 ' 1) . 
q i 2 " q 42 “ ,, “777 zrz 

(I'Y) (1-yC y) 

2 2 

11 _ 11 _ Y q 3 

q 13 ' " q 43 r n w-, c 2 . 

(1-y)(1-Y-£ y) 

nil = -nil = ^*3 

14 44 (1-y) (1-y-c 2 y) . 

qll f K 2 +1 ^ 

21 ~T~ 

a-C Y 

q 12 = " q 42 = 2 y/(1-Y-5 2 y) 
q 13 -q 43 -y/(1-Y-5 Y) 


q 21 = y (°l + IV (a ' 


.12 _ A ,, 2 . 
q 22 " (“ " 5 YJ 


*11 = -y(Oj * 1)/Ca - 5 2 f) 


q 21 = ‘ *V C “ ‘ ^ 


- 2 T /a-Y-E 2 T ) 


q 22 — y/(1-Y-? 2 y) 
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0 0 


q 2 2 \ • - a * K 2 )o/(« - s 2 v) 


<Jato 


q 2 J 0 0 

0 0 0 


0 


0 


qjl = (2a + K x - 1)/(1-Y-5 2 Y) 
q 43 = Y^/fl-Y-C^Y) 


0 





q 44 = yA 3 /(1-y-5 2 y) 


(6.16) 


In order to extract the subsystem matrices and p£ independent 
of the coupling parameter c and obtain the decomposition of (6.8) into 
(6.13-6.14) it is necessary to use the following identities: 


1 - 1 . g 2 Y 

a - 5 Y a (a - K y) 

L_ = -1-. + g 2 Y 

Y - 1 + 5 2 y y 1 (1-y) (y-1 + C 2 y) 


(6.17) 


The structural configuration of the system S as composed of the two 
subsystems and and the interconnections between them through the 
coupling parameter £ can then be depicted as in Fig. 2. 

It is obvious that the system of Fig. 2 becomes that of Fig. 3 when 

5 = 0 . - i ■ .. •- . • v ■ ■ ••• • - • 

On the basis of the physical characteristics of the Skylab given in * . 
Section 5, the matrices Q... (£) (1; j =1, 2) of (6.16) can be made indepen- 
dent of E, and denoted by . That is accomplished by neglecting the term 

£ 2 Y = 8.5 x 10' 4 (6.18) 
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5 f 0 



Fig. 6. 2 

in (6.16), with respect to the terms 

1 - y = 0.803 and a - 5.52 . 


Fig. 6. 3 


(6.19) 


After this simplification, the numbers (i = 1, 2) of the norm of 


the coupling matrices (L ^ can be computed using 


ly - [A(Qy Qy)] 1/2 . i, j = 1, 2 


( 6 . 20 ) 


and the aggregation matrix A = (a^) defined by (3.15), becomes a function 
of the coupling parameter £ only. 

The Computer Program given in the Appendix, is used to find: 
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T 1/2 i = 1 2 

a. Subsystem Liapunov functions = (xdho) ' ’ ’ ’ 

of (3.4); 

b. Numbers , i = 1, 2 ; j = 1, 2, 3, 4 , of (3.7) ; 

c. Numbers , i, j = 1, 2, of (6.20); 

d. Aggregation matrix A of (3.14) as a function of the coupling 
parameter £ ; and 

e. Solution of the stability inequalities (4.7) in terms of the 
maximum value £^ of £ . 

Subsystem Liapunov's functions v^ , i = 1, 2, are obtained by solving 
the Liapunov matrix equations 

pTh. + H.P = -G. , i = 1, 2 (6.21) 

using the direct method of solution as described in the Appendix. 

The choice of the 4x4 symmetric matrices = I i = 1, 2, where 
I is the 4x4 identity matrix, yields the positive definite 4x4 sym- 


metric matrices FL , i = 1, 2, 

and establishes global asymptotic stability 

of the subsystems , i = 1, 2 

• 


Then, to construct the aggregation matrix A , 

the following numbers are 

computed: 



XOy = 10.0811 , 

A (H^ » 5162.7539 

, A(G X ) = 1 

X(H 2 ) = 0.4176 

A(H 2 ) = 436.3635 

x(c 2 ) = ! 

n ll = ^*1750 

n 2 -|_ = 0.6462 


n^2 = '71.8523 

n 22 = 20.8893 
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n 13 = ^-0069 


n 14 = 1626.0258 


"23 

n 24 


0.0239 

675.2502 


K n = 0.36465 


5 12 = 0.77665 


5 21 = 1.84785 5^2 = 0.54915' 


( 6 . 22 ) 


Finally, the aggregation matrix A is obtained as 


A = 


-0.96 x 10 -4 + 186. 755 2 


392.985 


The stability inequalities 


-0.96 x 10" 4 .+ 186. 755 2 < 0 


1954.265 


-11.45 x 10 -4 + 573. 865 2 


(6.23) 


-0.96 x 10" 4 + 186. 755 2 


392.985 


1954.265 


-11.45 x 10' 4 + 573. 865 2 


> 0 . 


' (6.24) 


are satisfied for all 5 such that 


0 < 5 i 5 m = 0.38 x 10 


-6 


(6.25) 


Hie obtained range (6.25) of the coupling parameter 5 is small due to 
the conservativeness of the stability procedure. However, the estimate 5 m 
could be considerably increased by a proper choice of the matrices , 
i = 1, 2, in (6.21). A meaningful optimization problem can be formulated as 
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the maximization of £ over all matrices G. . Future effort should be 

m 1 

directed toward a solution of this optimization problem, which can provide 
important information about the trade-off that exists between the degrees 
of the subsystem and the over-all system stability. 


7. ACTIVE CONTROL 


A mission requirement of the spinning Skylab is that the 3-axis be pointed 
at the sun. In order to inert ially fix the 3 -axis in presence of disturbance 
torques, attitude control torques must be applied to the vehicle [1, 12]. The 
control torques depend on error signals that are proportional to the angle be- 
tween the 3-axis and the solar vector. Sun sensors and rate gyros on the pre- 
sent Skylab can readily provide the control signals 4>p <J> 2 > and w 2 shown 
on Fig. 5.2. 

The linear control is postulated [12] as 


T = ot<i> + 


(7.1) 


T T 

where T = [T^ T 2 T^] is the vector of control torques; $ = [<f>^ <p^ i^] is 

T 

the vector of angular rotations; uj = [w^ w^ w^ + 0] is the vector of angular 
velocities; a , 0 are 3><3 matrices 


a ll 

a 12 

1 

o 


*11 

*12 

o 

a 21 

a 22 

0 

, B = 

*21 

*22 

0 

O 

0 

0_ 


0 

0 

*33 


and kinematic relationships are 



“1 

0 

0“ 


~0 

a 

0~ 

a) - 

0 

-1 

0 

+ 

-n 

0 

0 


_0 

0 

-1_ 


_0 

0 

0_ 


(7.2) 


r 


(7.3) 
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The control law in this study is chosen as 


a 12 I l fi ’ 


all other ou ^ = 0 


= 1^6 , = - I^fip, all other = 0 


(7.4) 


so that the normalized control torques v = 
Tj/I n 2 -] T are 


1 V 1 v 2 V 3 1 T * 'V 1 / W 2 


V 1 = (e + 6)4 2 - 6 (j> 1 


Vi = 0 


v, = p 4- 


(7.5) 


Referring to equations (5.3) and (5.4) , the control torque is used to sta- 
bilize the subsystem (wobble motion), and the torque is used to sta- 
bilize the subsystem S 2 (spin motion) . 

In (7.5), e, 6, p are control parameters to be selected in the stabili- 
zation process. 

Upon introducing these transformations the linearized equations of motion 


become : 


wobble 

motion 




fl 


I 


M 


t fl 


1 - (1 + K 1 )4>2 " ^1^1 " Y ^ y 3 +y 3^ + SY^y-j+y^y^ + ( e+6 ) < < > 
(1-*-K- l ) 4> 1 + f> 2 + ' CY(w 1 -y 1 -2y 2 ) = 0 


2 



= 0 


'h ■ *1 + y 3 + A 3 y 3 + (a 3 + 1)y 3 = 0 


(7.6) 
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spin 

motion 


B4>-r + yOh - 2y~) + p<}>, = 0 


ft ft ff 


-2e,4>i - C4>2 + $ t, + y l + A l y l + °l y l 


2y 2 + “ 0 


£<!>! - 2 C 4>2 ' + 2 4> 3 + 2 u ^ + u 2 + A 2 y 2 + ^°2 * 1 ^ y 2 = 0 


(7.7) 


The state space representation of the overall system S described by 
(7.6-7), is obtained using the same method as in the Passive Control case and 
is : 

S: x (x) = P x(t) , (7.8) 

where the state 11-vector x(t) is chosen as 

» » ? i » * •p 

X ( T ) — ($j» < f > 2> ^2* y 3> < f > 3» y ]_» ^ 2 * y 2^ (7*9) 

and the 11 * 11 matrix P is given in Fig. 7.1. 

The system S of equation (7.8) can be decomposed into two interconnected 
subsystems described by: 

(7.10) 

(7.11) 


motion S l : X 1 ^ = P 1 X 1^ T ^ + + 

spin S 2 : x 2 (t) = P 2 x 2 (t) + 5Q 21 C5JX 1 (t) + l Q 2 2(O x 2( t ) 
motion 
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using the same procedure outlined in the previous section. 

The state vectors x(t), (t), (t) of the system S and two subsys- 

tems and are 


x(t) = 


’x n (t)‘ 


x, (t) 



(7.12) 


In (7.10-11) the 6^6 and 5x5 matrices P^, P 2 correspond to the sub- 
systems and and 6 x 6, 6 x 5, 5 x 6 and 5 x 5 matrices Q-jjGi), 
Q 12 (?), Q 21 (C), Q 22 (5) represent the interconnections between the two subsys- 
tems: 


0 0 0 1 0 0 

0 0 0 0 1 0 

0 0 0 0 0 1 

1 1 1 1 1 1 

P4I P42 p 43 p 44 p 45 P 46 

0 p 52 0 p 54 0 0 

P 61 P 62 P 63 P 64 P 65 P 66 


P 41 = Ck 1 +y)/(i-y) 

p 42 “'P62-"' _(e+6)/(1 " Y) 

p 43 = " ^3/a-Y) 

P44 = Pg4 = «/d-Y) 

p 45 = P 65 = Cl+y/Cl-Y) 
p 46 = P 66 = -YA 3 /d-Y) 


P 52 = “ K 2 

Pm’V 1 

Pei = ( 1 +k i)/Ci-t) 

P63 '- f 4 t i-D/Ci-r) 
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2 2 2 
P 11 P 12 P 13 

0 

0 




0 

0 1 

0 

0 



P 2 

2 2 2 

P 31 P 32 P 33 

0 

2 




0 

0 0 

0 

1 




-2 

0 -2 

P54 

2 

P55 



2 

Pll 

= . 

■ p/(6-y) 

4 = 

- Ba^/(6-Y) 


2 

Pl2 

? Y o^/(3-y) 

2 

p 33 

- ba^cb-y) 


2 

Pl3 

= YA 1 /(e-Y) 


- Cp 2 -1) 


2 

P 3 1 

= p/(3*y) 

P 55 * 

"*2 , 





~ 0 0 

0 

0 

0 

0 ' 



0 0 

0 

0 

0 

0 



0 0 

0 

0. 

0 

0 

Q H(0 

= 

11 11 

q 41 q 42 

q 43 

4 

4 

4 



n 11 

0 q 52 

0 

4 

0 

0 



11 11 
q 61 %2 

<& 

4 

4 

q 66 

4 

11 YQC-J+2Y-1) 

q 61 n . n 

(1-Y) (1-Y" 

?7) 


4= 

%S = 

4 

= qjl = - X<fr«L-. 

(1-y) (1-y 

-£ 2 y) 


4 = 

^6 = 


(K^y-iOy 
(1-y) (1-Y-C^y) 

(1-y) (1-y-C Z y) 
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11 - „11 — 

q 43 “ q 63 


y 


(l-y) (1-yC r) 


11 _ „H _ Y5 

q 44 “ q 54 " 


(1-y) (l-y-rr) 


11 

q S2 


q 54 


y(y i) 

a(l-y 3 )-C Z Y 

Y(3-K 2 ) 

a(l-Y 3 )'C^Y 


Qi 2 (5) 


0 

0 

0 

12 

q 41 

12 

q 51 

12 

q 61 


0 

0 

0 


12 

q 52 


0 

0 

0 

0 

3 12 

^53 


0 

0 

0 

12 

q 44 


12 

V 


0 

0 

s 

12 

q 45 


12 

q 65 


12 _ „12 _ 2y 
q 41 “ q 61 


T 


1-Y-C Y 

12 _ .12 _ Ya 2 
q 44 q 64 " ^7^7 


1.2 _ y 3 p 
q Sl " 


a(l-Y 3 )-£ 'Y 

.12 _ Y( l-Y 3 +a l) 
52 a(l-y 3 )-5 Y 


12 _ 12 
q 45 " q 65 


yA 2 

i-y-c 2 y 


12 

q 53 




a(l-Y 3 )-C Y 


22 

q ll 


? 22 

*12 


22 

hz 


Q 22 CO 


22 

q 31 


22 

q 32 


j 22 

*33 


22 

^51 


0 

0 


0 

0 

51 


0 

22 

q S4 


22 

q 55 



22 

<*11 


22 

^12 


22 

*13 


PY, 


(l'Yj) [“(l'Yj)'C y] 

YYgCo^-Yg) 

(1-Y 3 )[a(l-Y 3 )-C 2 Y] 

^3 A 1 

(1-Y 3 )[a(l-Y 3 )-C 2 Y] 


22 

q 33 


,22 _ 


*51 


22 

*S4 


YA 1 

(1-Y 3 )[a(l-Y 3 )-C 2 Y] 

2y _ 

1-Y-C 2 Y 

. 2 
Y°2 

(1-Y-C 2 Y) 


22 

q 31 


PY 


ftd-Y 3 )[ad-Y 3 )-E Y] 


22 

*55 


ya, 


1-y-5 2 y 


22 Y(a 1+ 1-Y 3 ) 

q 32 


Q 21 (C) " 


d-Y 3 )[a(l, 

-r 3 )-c 2 rl 



1 

o 

21 

q 12 

0, 

„21 

q 14 

0 

0 

0 

0 

0 

„21 

q 32 

0 

21 

q 34 

0 

0 

0 

0 

21 - . 
q 5i - 

*21 

q 52 

21 

q 53 

21 

q 54 


0 

0 

0 

0 


21 

q 55 


0 

0 

0 

0 

i 21 

*56 


21 Y 3 a(K 2 + 

q 12 


(1-y 3 )-C y 
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1-Y-C Y 


21 _ 

q 5i 7 c 2 
1-Y-C Y 


21 _ yA 3 


<56 : ~r 

l-Y-C Y 


(7.13) 
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The following identity relationships were used in order to get the sub- 
system matrices and P£ independent of the coupling parameter £ : 


a(l-Y 3 )-5 2 Y " a(1 " Y 3 ) a(l- Y 3)[aCl-Y 3 )-5 2 Y] 

1 JL + c 2 y 

i-y-5 2 y 1 Y (1 - y)(1-y-5 y) 


(7.14) 


The graphical interpretation of the interconnected subsystems and 

S 2 is the same as in Fig. 2 and 3 of the previous section. 

Again, on the basis of the physical characteristic of the Skylab given in 

Section 5, the matrices Q. . (O (i, j = 1, 2) of (7.13) can be made independent 

2 

of 5 and denoted by Q. . . This is obtained by neglecting the term 5 y = 

-4 

8.5 x 10 with respect to the terms 


1-Y = 0.803 and cx(1-y 3 ) = 5.33 


(7.15) 


Using the following specific values of the control parameters e, 6, p: 


e = 2.0 5 = - 1.0 p = 1.0 


(7.16) 


the same computational algorithm as in the Passive Control case is applied and 
the computer results are shown in the Appendix. 

The choice of the 6^6 and 5 x 5 identity matrices for the G matri- 
ces of the first and second subsystem results in 6 x 6 and 5x5 positive 
definite matrices H^, FL, and establishes the global asymptotic stability of 
the subsystems. 

In order to construct the aggregation matrix A , the following numbers 
are computed: 
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XCHj) = 0.4947 , A^) = 44.2273 , A(G X ) = 1 


X(H 2 ) = 0.4176 , A(H 2 ) = 427.5991 , X(G 2 ) = 1 


= 0.7034 
n^2 = 6.6503 
= 0.0751 
n 14 = 62.8749 

C n = 0.7897C 2 
C 21 = 2.7088C 


n 2 i = 0.6462 
n 22 = 20.6784 
n 2 2 = 0.0241 
n 24 = 661.6879 

C 12 = 314.5269c 
C 22 = 222.4041C 2 


(7.17) 


The aggregation matrix A is obtained as 


A = 


-11.29 x 10“ 3 + 70.59C 2 


2543. 19C 


30603. 29C 


-1.16 x 10' 3 + 227734. 68C 2 


(7.18) 


The stability inequalities 


-11.29 x 10" 3 + 70.59C 2 < 0 


-11.29 x 10 _3 + 70.59C 2 
2543. 19C 

are satisfied for all C such that 
0 <_ C < C m = 0.41 x lo" 6 . 


30603. 29C 


-1.16 x 10‘ 3 + 227734. 68C 2 


> 0 

(7.19) 


(7.20) 
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The obtained interval (7.20) of the coupling parameter £ for which the 
overall system is globally exponentially stable is relatively small due to the 
following reasons : 

1. The inherent conservativeness of the stability analysis; 

2. The choice of the matrices , i = 1, 2, is not the 

"best" regarding the value of £ ; and 

m 

3. The freedom in the choice of the control parameters e, 

6, p is not used to the full extent. 

Future research should elaborate on the points 2 and 3 and formulate an 
optimization problem: Maximization of £ m over the elements of the matrices 

G^ , i = 1, 2, and the control parameters e, 6, p . Nonsystematic numerical 
experimentation with parameter values indicated a possibility of considerable 
improvements in the value of £ m . 
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CONCLUSION 


A decomposition- aggregation method is outlined for stability analysis of 
large-scale dynamic systems. The method takes advantage of special structural 
features of the complex systems to reduce the memory and computational time re- 
quirements when the stability analysis is carried out by machine calculations. 

By utilizing the system structure in the decomposition procedure, the decompo- 
sition-aggregation method makes explicit important structural properties of 
the system. Furtheimore, the method is suitable for accommodation of nonlin- 
earities either in the subsystems or in their interconnections. However, the 
method is inherently conservative since a series of approximations are involved 
in establishing the sufficient conditions for stability. Therefore the success 
of the method should be judged satisfactory to the extent that the conservative- 
ness of the results is outweighed by the computability of the method and the in- 
sight that the method provides into the structural properties of complex dynamic 
systems . 

The decomposition-aggregation method is applied to the dynamic model of a 
spinning Skylab. After the model is decomposed into the wobble and spin subsys- 
tems, both the passive and the active control are considered. Such decomposi- 
tion made an important structural parameter to appear as an interconnection 
parameter of the two subsystems. Subsequent stability analysis was aimed at 
estimating the interval of the parameter for global exponential stability. The 
obtained estimates turned out to be relatively conservative since the flexibil- 
ity of the decomposition aggregation method was not used to the full extent. 
Moreover , several physical constraints of the control should have been removed 
in order to increase the degree of stability of the subsystems and achieve a 
higher degree of stability on the over-all system level. Since the outlined 
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decomposition-aggregation analysis is completely computerized, the proposed 
improvements can be readily incorporated in the present analysis scheme to 
yield a flexible and powerful method for stability analysis of large-scale 
linear and nonlinear dynamic systems. 
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APPENDIX 


COMPUTER PROGRAMS 



DESCRIPTION OF THE COMPUTER PROGRAM 

The program for analyzing the passive stability of the Skylab using de- 
composition-aggregation method and vector Liapunov functions is realized on 
the IBM 1130 computer (16 K memory) in FORTRAN language. 

In the following description the subroutines: LOC, MSTR, EIGEN, SIM) 
and MATA are IBM supplied subroutines (from IBM 1130 Scientific Subroutine 
Package). Storage compression feature was used for handling the arrays in 
these subroutines. The three modes of storage are termed general, symmetric , 
and diagonal. General mode is one in which all elements of the matrix are 
in storage. Symmetric mode is one in which only the upper triangular portion 
of the matrix is retained columnwise in sequential locations in storage. 

(The assumption is made that the corresponding elements in the lower triangle 
have the same value) . Diagonal mode is one in which only the diagonal ele- 
ments of the matrix are retained in sequential locations in storage. (The 
off diagonal elements are assumed to be zero) . This capability has been im- 
plemented using vector storage approach. 

The names of the variables in the mainline program were chosen such that 
they either completely match to the notation throughout the text, or in cases 

where it was impossible to strongly indicate what was meant (e.g. , tu-, and 

0 

ETA 1(1)). 

The flowchart, description of the subroutines and the computer programs 
are included. 
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Write the matrix 
L G / 
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Subroutine LOC 


Purpose: Compute a vector subscript for an element in a matrix of specified 

storage mode. 

U*age: CALL LOC (I, J, IR, N, M, MS). 

VescAlptlon o & pasumetens : 

I - Row number of element 
J - Column number of element 
IR - Resultant vector subscript 
N - Number of rows in matrix 
M - Number of columns in matrix 
MS - One digit number for storage mode of matrix: 

0 - general , 1 - symmetric , 2 - diagonal 

Method. : 

MS = 0 Subscript is computed for a matrix with N*M elements in storage 
(general matrix) 

MS = 1 Subscript is computed for a matrix with N*(N+l)/2 in storage 
(upper triangle of symmetric matrix) . If element is in lower 
triangular portion, subscript is corresponding element in upper 
triangle. 1 - - : - 

MS = 2 Subscript is computed for a matrix with N elements in storage 
-(diagonal elements of diagonal matrix). 'If element is not. on 
diagonal (and therefore not in storage) IR -is set to zero. 

Subroutine MSTR 

Purpose: Change storage mode of a matrix 

ll&age: CALL MSTR (A, R, N, MSA, MSR) 

Description oj parameters : . ..... 

A - Name of input matrix 
R - Name of output matrix 

- ■ i . ■ '■ , "ti 

N - Number of rows and columns in A and R 
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MSA - One digit number for storage mode of matrix A 
0 - general , 1 - symmetric , 2 - diagonal 

Rmcuiki : Matrix R cannot be in the same location as matrix A . Matrix 

A must be a square matrix. 

Subiootine and function iubpAogAam A&quuAed : LOC 

M eXhod : Matrix A is restructured to form matrix R . 

MSA MSR 
0 0 

0 1 

0 2 

*\ 

1 0 

1 1 

. 1 2 

2 0 

' • ■-2 . . - '1 

2 2 

SubAoutim LOCI 

The same as subroutine LOC except that the symmetric mode is to be con- 
sidered one in which only the upper triangular portion of the matrix is- re- 
tained but aom-uxc6 <l in sequential locations in storage. • 

SubAoutiw MSTRl | 

The same as subroutine MSTR with the same remakr about symmetric mode 


Matrix A is moved to matrix R 

The upper triangle elements of a general matrix are used to 
form a symmetric matrix 

The diagonal elements of a general matrix are used to form 
a diagonal matrix 

A symmetric matrix is expanded to form a general matrix 
Matrix A is moved to matrix R 

The diagonal elements of a symmetric matrix are used to foim 
a diagonal matrix 

A diagonal matrix is expanded by inserting missing zero ele- 
ments to form a general matrix 

-A diagonal matrix is expanded by inserting missing zero ele- 
ments to form a symmetric matrix 

Matrix A is moved to matrix R 
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Subroutine EIGEN 


Purpos e : Compute eigenvalues and eigenvectors of a real symmetric matrix 

Usage: CALL EIGEN (A, R, N, MV) 

Description oh parmeten> : 

A - Original matrix (symmetric) destroyed in computation. Resultant 
eigenvalues are developed in diagonal of matrix A in descending 
order 

R - Resultant matrix of eigenvectors (stored columnwise, in same se- 
quence as eigenvalues) 

N - Order of matrices A and R 

MV - Input code: 

0 - Compute eigenvalues and eigenvectors 

1 - Compute eigenvalues only (R need not be dimensioned but must still 

appear in calling sequence) 

Remarks : Original matrix A must be real symmetric (storage mode = 1) . Matrix 
A cannot be in the same location as matrix R . 

Method : Diagonal ization method originated by Jacobi and adapted by von Neumann 

for large computers. 


Subroutine SIMQ 

Purpose: Obtain solution of a set of simultaneous linear equations, AX = B. 

Usage: CALL SIMQ (A, B, N, KS) 

Description ofi parameters : 

A - Matrix of coefficients stored solumnwise. These are destroyed in 
the computation. Matrix is N by N . 

B - Vector of original constants (length N). These are replaced by 
final solution values, vector X . 

N - Number of equations and variables. N must be greater than 1. 

KS - Output digit: 

0 - For a normal solution 

1 - For a singular set of equations 
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Remarks : Matrix A must be general. If matrix is singular, solution values 

are meaningless. 

Method : Method of solution is by elimination using largest pivotal divisor. 


Subroutine MATA 

Purpose.: Premultiply a matrix by its transpose to form a symmetric matrix. 

U&age: CALL MATA (A, R, N, M, MS) 

Description oj parameters •• 

A - Name of input matrix 
R - Name of output matrix 
N - Number of rows in A 

M - Number of columns in A . Also number of rows and number of columns 
of R 

MS - One digit number for storage mode of matrix A 
0 - general , 1 - symmetric , 2 - diagonal 

Remarks Matrix R cannot be in the same location as matrix A . Matrix R 
is always a symmetric matrix with a storage mode = 1. 

Subroutine and function subprogram required: LOC 

Method : Calculation of (A transpose A) results in a symmetric matrix re- 

gardless of the storage mode of the input matrix. The elements 
of matrix A are not changed. 

Subroutine LIAPU 

T 

Purpose: Expanding the Liapunov matrix equation P H + HP = -G into a system 

of linear algebraic equations. 

Usage: CALL LIAPU (P, L, N, M, V) 

Description oj parameters : 

P - Name of input matrix 

L - Auxiliary matrix consisting of integers and used to construct the 
matrix V of the corresponding system of linear equations 
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N - Order of matrices P and L 

M - Order of the system of linear equations, that is, order of matrix 
V . Also equal to N(N+l)/2 

V - Name of output matrix of transformed system of linear equations 

Rma/iki. : The transformed system of M linear equations is V-HR = -GR , 

where: 

GR - Vector formed from symmetric matrix G using the upper triangular 
elements of the matrix G stored Kow-mAZ 

HR - Vector of unknowns formed from unknown symmetric matrix H using 
the upper triangular elements of H stored Aow-uiciz 

Method : Expansion of the matrix equation P^H + HP = -G into N*(N+l)/2 

simultaneous equations, using the method developed in [ 7 ] and re- 
fined in [ 8 ] 
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04/24/73 


JOB NO. 023608 

// JOB T _ 

LOG DRIVE CART SPEC CART AVAIL PHY DRIVE 
0000 0001 0001 0000 

V2 MIO ACTUAL L6K CONFIG I6K 

// * 

// * PASSIVE STABILITY OF THE SPINNING SKYLAB ' 

// *• 

// FOR ' '• V 

♦ONE WORD INTEGERS ‘ 

♦LIST SOURCE PROGRAM t 

SUBROUTINE LOCI < I , J, IR,N,M,MS ) 

I X*I 
JX = J 

IF (MS-n 10,20,30 
10 IRX*N*( JX-D + IX 
GO. TO 36 

20 IF(IX-JX) 22,24,24 
22 I RX= J+ ( I - I ) ♦ ( 2*N- I ) / 2 
GO TO 36 

24 I RX= I ♦ ( J-I ) * ( 2*N-J ) ./ 2 
GO TO 36 
30 I RX=0 

IF(IX-JX) 36,32,36 
32 I RX= I X 
36 I R=I RX 
RETURN 
END 

FEATURES SUPPORTED 
ONE WORD INTEGERS 

CORE REQUIREMENTS FOR LOCI 
COMMON 0 VARIABLES 6 PROGRAM 130 

RELATIVE ENTRY POINT ADDRESS IS 0009 (HEX) 

END OF COMPILATION 

// DUP 

♦STORE WS UA LOCI 

CART ID 0001 DB ADDR 501F DB CNT 0009 
// EJECT 
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//FOR . _ _ 

♦ONE WORD INTEGERS 
♦LIST SOURCE PROGRAM 

SUBROUTINE MSTR1 ( A , R , N , MSA , MSR ) 
DIMENSION A(1),R(1) 

DO 20 I = 1 » N 
DO 20 J = 1,N 
IF (MSR) 5,10,5 
5 IF(I-J) 10,10,20 
10 CALL LOCK I , J , I R , N , N , MSR ) 

IF(IR) 20,20,15 
15 R ( I R ) =0 • 0 

CALL LOCK I , J, IA,N,N,MSA) 

IF(IA) 20,20,18 
18 RCIR)=A(IA) 

20 CONTINUE 
RETURN 
END 

FEATURES SUP PORT fD 
ONE WORD INTEGERS 

CORE REQUIREMENTS FOR MSTR1 
COMMON 0 VARIABLES 6 PROGRAM 

RELATIVE ENTRY POINT ADDRESS IS 0009 (HEX) 

END OF COMPILATION 

// DUP 

♦STORE WS UA MSTR1 

CART ID 0001 DB ADDR 5028 DB CNT 0008 
// EJECT 


110 
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// FOR - _____ 

* LIST SOURCE PROGRAM 
♦ONE WORD INTEGERS 

SUBROUTINE L I APU ( P , L , N, M , V ) 

DIMENSION L(4f4) ,P(4,4),V(10,10) 

M=N* ( N+ 1 ) / 2 
K*0 

DO 10 1=1, N 
DO 10 J=I,N 
K = K + 1 
L ( I , J ) =K 

10 L ( J » I ) =K 

DO 11 1 = 1, M 
DO 11 J= 1 » M 

11 V(I,J)=0. 

DO 12 1 = 1, N 
DO 12 J=1,N 
DO 12 K= 1 , N 
1 1 =L ( I » K ) 

I J = L ( J » K ) 

12 VIII ,IJ)=P(J,I)+V(II,IJ) 

DO 13 I = 1 , N 

DO 13 J= 1 , M 
I L=L ( I , I ) 

13 V( IL , J ) =2 « *V ( I L, J ) _ 

RETURN 

END 

FEATURES SUPPORTED 
ONE WORD INTEGERS 

CORE REQUIREMENTS FOR LIAPU 
COMMON 0 VARIABLES 10 PROGRAM 

RELATIVE ENTRY POINT ADDRESS IS 0011 (HEX) 

END OF COMPILATION 

// DUP 

♦STORE WS UA LIAPU 

CART ID 0001 DB ADDR 5030 DB CNT 0012 


280 


// EJECT 
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// FOR . 

*♦ PASSIVE STABILITY OF THE SPINNING SKYLAB 
♦ONE WORD INTEGERS 

♦LIST SOURCE PROGRAM • - 

♦ I OC S ( CARDj 1403 PRINTER 
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04/24/73 PASSIVE STABILITY OF THE SPINNING SKYLAB 

__ REAL II, 12,13, H ASS, K1,K2 

DIMENSION L!4,4) ,P{4,4) , Q ( 4 , 4 ) , G ( 4 , 4 > , H< 4 , 4 ) , GC( 1 0 ) , HC ! 10 ) ,GR 110 ) 
DIMENSION V! 10, 10) ,GD<4) ,HD<4) ,E(4) 

DIMENSION ETA1!2),ETA2!2),ETA3!2),ETA4!2),A!2,2) 

EQUIVALENCE <P<l,l) f Q(l,n f G<l,l>,H<l,in,(GC(l>,HC(l>> 

100 FORMAT (12) 

101 FORMAT ( / ) 

102 FORMAT!//) 

103 FORMAT! 10F12. 4) 

C READ N - THE ORDER OF THE SUBSYSTEM 

RE AD ( 2 » 100 ) N 

C READ THE PHYSICAL PARAMETARS OF THE SKYLAB 

RE AD ! 2 ,50 ) 1 1 , 12 , 1 3 , G2 , MASS , EK1 , EK2, OMEGA 
50 FORMAT ! 8F 10. 0 ) 

C COMPUTE THE NORMALIZED PARAMETARS 

K 1 = ( I 2- 1 3 ) / 1 1 
K2= ( 1 3-1 1 ) / I 2 
ALPHA=(1.0+Kl)/tl.0-K2) 

GAMA=(2«*MASS*G2*G2) /I1 
SI GS l = EKl/ ( MASS*OMEGA*OMEGA ) “ 

S I GS2=EK2/( MASS* OMEGA* OMEGA) 

SIGM1=SQRT ( SIGS1 ) 

SIGM2=SQRT!SIGS2) 

DEL1=0.04*SIGM1 

DEL2=0.04*SIGM2 

GAMA 1 = I . O-GAMA 
GAMA2=GAMA1*GAMA1 
WRITE!5,301) 

301 FORMAT! *1' ,10X, 'STABILITY ANALYSIS OF LARGE SCALE SYSTEM USING DEC 
lOMPQSI T I ON METHOD AND LIAPUNOV FUNCTIONS’,//) 

S U B S Y S T E M A N A _ L Y S I S 

DO 9 K = 1 , 2 
DO 5 1=1, N 
DO 5 J = 1 , N 

5 P(I,Ji=0.0 

IF(K-l) 7,6,7 

6 WR I TE ( 5 , 302 ) 

302 FORMAT!' THE FIRST SUBSYSTEM MATRIX P IS',/) 

C COMPUTE THE ELEMENTS OF THE FIRST SUBSYSTEM MATRIX PI 

P( 1,2)=!K1+GAMA)/GAMA1 

P< lj3) = (GAMA*SIGSl)/GAMAl_ _ 

P( 1,4) =!GAMA*DEL1) /GAMA 1 

P(2,1)=K2 

P(3,4)=1.0 

P!4,2)=-!K1+1.0)/GAMA1 
P!4,3)=(GAMA-SIGS1-1.0) /GAMA1 
P!4,4)=-DEL1/GAMA1 
GO TO 8 

7 WR I TE ( 5 , 303 ) 

303 F0RMATC1THE SECOND SUBSYSTEM MATRIX P IS’,/) 

C COMPUTE THE ELEMENTS OF THE SECOND SUBSYSTEM MATRIX P2 

P! 1,2)=1.0 
P ( 2 , 1 ) =-S I GS 1 
P < 2 , 2 ) =-DEL 1 
P ( 2 , 4 ) =2 • 0 
P<3,4)=1.0 
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P(4, 2 )=-2. 0 ■ i 

P ( 4 » 3 ) * 1 .O-SI GS2 
...P (4 f 4 ) =-0EL2 

C ' WRITE THE CORRESPONDING SUBSYSTEM MATRIX 
8 DO 16 1 = 1, N 

16 WR I TE ! 5 * 103 ) ( P ( I , J ) , J = 1 , N ) 

WR1T£15jlIQJJ. 

FROM THE MATRIX P AND THE MATRIX EQUATION PTRANSPOSE+H + H*P = -G 
COMPUTE THE MATRIX V OF THE TRANSFORMED SYSTEM OF LINEAR EQUATIONS 
CALL LIAPU(P,L,N,M,V) 

WR I TE ( 5 , 304 ) 

304 FORMAT! • THE MATRIX V OF THE CORRESPONDING SYSTEM OF LINEAR EUUATI 

IONS ISM 

WRITE (5,101) 

WRITE V MATRIX 
DO 14 1=1, M 

14 WRI TE ( 5 , 103 ) ( V ( I , J ) , J= 1 , M ) 

WRITE(5,101) 

READ £R - VECTOR CONSISTING OF UPPER TR I ANGULAR ELEMENTS 
OF MATRIX G STORED ROWISE 
READ ( 2 , 104 ) GR 

104 FORMAT! 16F5.0) 

OBTAIN THE GENERAL MATRIX G (STORAGE MODE 0) 

CALL MSTR1 (GR,G,N,1 ,0) 

WRITE(.5 l 109) 

109 FORMAT! • THE POSITIVE DEFINITE SYMMETRIC MATRIX G IS',/) 

DO 19 1=1, N 

1 9 WR ITE!5,103) l G ! I , J ) , J= 1 , N ) 

OBTAIN GC - VECTOR CONSISTING OF THE UPPER TRIANGULAR ELEMENTS 
OF THE MATRIX G STORED COLUMNWISE 
CALL MSTR!G,GC,N,0,1) 

SOLVE FOR THE EIGENVALUES OF THE MATRIX G 
CALL E I GEN ( GC , 0 , N , 1 ) 

STORE THE EIGENVALUES OF THE MATRIX G IN THE VECTOR GD 
CALL MSTR!GC,GD,N,1,2) 

WR I T E ! 5 , 10 1 ) 

WRITE (5 ,201 ) 

201 FORMAT! » THE EIGENVALUES OF THE SYMMETRIC MATRIX G ARRANGED IN DEC 
1 RE AS I NG ORDER ARE* ,/) 

WRITE THE EIGENVALUES OF THE MATRIX G 
WRITE!5,103) GD 
WRITE!5,102) 

SOLVE THE SYSTEM OF LINEAR EQUATIONS V*HR=GR, RESULT IS IN GR 
DO 17 1=1, M 

17 GR! I ) =-GR ! I ) 

CALL SIMQ! V,GR,M,KS) 

IF(KS-l) 30,20,30 

20 WRITE!5,105) 

105 FORMAT! • SINGULAR CASE') 

GO TO 15 

FORM GENERAL MATRIX H FROM GR 
30 CALL MSTR1!GR,H,N,1,0> 

WRITE!5, 110) 

110 FORMAT!' THE LIAPUNOV MATRIX H FOR THE SUBSYSTEM IS',/) 

DO 21 1=1, N 

21 WR I TE ! 5, 103 ) < H ! I , J ) , J= 1 , N ) 

OBTAIN HC - VECTUR CONSISTINGUF THE UPPER TRIANGULAR ELEMENTS 
OF THE MATRIX H STORED COLUMNWISE 
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CALL M STFUH.HC.N.O, I) ■ •_ 1 

C SOLVE FOR THE EIGENVALUES OF THE MATRIX H 

CALL EIGEN(HC,D,N,I) 

C STORE THE EIGENVALUES OF H IN HD 

CALL MSTR(HC,HD,N,1,2) 

WRITE (5, 101) 

WRJTTE ( 5j 202 ) _ 

202 FORMAT! • THE EIGENVALUES OF THE SYMMETRIC MATRIX H ARRANGED IN DEC 
1REASING ORDER ARE',/) 

C WRITE THE EIGENVALUES OF THE MATRIX H 

WR I TE ( 5 » 103 ) HD 
WR I TE ( 5 » 101 ) 

C CHECK IF ALL EIGENV ALUES ARE PO SITIVE _ 

DO 22 1=1, N 
I F ( HD ( I ) ) 23,23,24 

22 CONTINUE 

23 WR I TE ( 5 , 203 ) 

203 FORMAT!' THE SUBSYSTEM" I S NOT ASYMPTOT ICALLY STABLE SINCE H MATRIX 

1 IS NOT POS ITIVE DEFINIT E* ) _ _ . 

GO TO 15 ““ . 

24 WR I TE ! 5 , 204 ) 

204 FORMAT! • ALL THE EIGENVALUES OF H MATRIX ARE POSITIVE AND THE SUBS 
1YSTEM IS ASYMPTOTICALLY STABLE') 

WRITE (5, 101) 

c _ COMP U TE Th e FOUR E STIMATES FO R THE SUBSYSJEM 

ETA1 ( K ) =SQRT ( HD ( N ) ) 

ETA2!K)=SQRT <HD( 1) ) 

ETA3 ( K ) =0. 5*GD( N ) /ETA2 ( K ) 

ETA4!K)=HD( 1 J/ETAUK) 

WR I TE ( 5 , 205 ) 

205 FORMAT! ' THE FOUR ES TIMATES FO R THE SUBSYSTEM ARE 'j / ) ; 

WRITE! 5, 103) E T A1 ( K ) , ETA2 ( K ) , ET A3 ( K ) , E T A4 i K ) 

9 CONTINUE 

FINDING THE NORMS OF THE COUPLING MATRICES 

WR I T E ( 5 , 2 06 ) - ■_ 

206 FORMAT! • 1 ' ,10X, 'ESTIMATING THE NORMS OF THE INTERCONNECTING MATRIC 
1 ES ' ) 

WR I TE ( 5 , 1 02 ) 

DO 52 K= 1 , 4 
DO 42 1 = 1 ,N 
DO 4 2 J = 1 , N 

42 Q ! I , J ) =0. 0 
IF(K-l) 44,43,44 

43 WR ITE ( 5 , 207 ) 

207 FORMAT !• THE SELFCOUPLING MATRIX Qll IS',/) 

C COMPUTE THE SELF-COUPLING MATRIX Qll 

Q< 1,2)=GAMA* (2.0*GAMA+K1-1.0)/GAMA2 
Q ( 1 , 3 ) =S IGS1*GAMA*GAMA/GAMA2 
0! 1,4)=(DEL1*GAMA*GAMAJ/GAMA2 
Q!2, 1) =GAMA*( K2+ 1.0) /ALPHA 
Q ( 4 , 2 ) =-Q (1,2) 

Q ( 4 , 3 ) =-Q (1,3) 

Q(4,4)=-Q( 1,4) 

GO TO 49 

44 I F ( K-2 ) 46,45,46 

45 WR I TE ( 5 ,208 ) 
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_ 208 FORMAT ( ///// , ' THE INTERCQUPLING MA TRIX Q12 IS*./) _ 

C COMPUTE THE INTER-COUPLING MATRIX Q12 
Q( 1,2)=2.0*GAMA/GAMA1 
Q< 1,3)=-GAMA/GAMA1 
0(2, 1)=GAMA*(S IGSl + l.O) /ALPHA 
Q(2»2)=GAMA*DEL1/ALPHA 

Q ( 4 .2)=-Q( 1.2) _ _ 

Q (4* 3 ) =-Q ( l » 3 ) 

GO TO 49 

46 IF ( K — 3 ) 48,47,48 

47 WR I T£ ( 5 » 2 10 ) 

210 FORMAT { • 1 THE SELFCOUPLING MATRIX Q22 IS',/) 

C COMPUTE THE SEL F-COUPL ING MATR IX Q22 __ 

Q ( 2 , 1 ) =-GAMA* ( S IGS1+ 1 . 0 ) /ALPHA 
Q(2,2)=-GAMA*DEL1/ALPHA 
Q ( 4 , 2 ) =2 • 0*G AMA/GAMA1 
0(4,3)=-GAMA/GAMAl 
GO TO 49 

48 WRIT E (5,21 1 ) _ ' 

211 FORMAT (/////, • THE INTERCQUPLING MATRIX Q21 IS',/) 

C COMPUTE THE INTER-COUPLING MATRIX Q21 

Q ( 2 , 1 ) =- ( 1 . 0+K2 ) 

Q(4,2)=(2.0*GAMA+K1-1.0)/GAMA1 

Q(4,3)=GAMA*SIGS1/GAMA1 

0(4,4)=GAMA*DEL1/GAMA1 

C WRITE THE CORRESPONDING COUPLING MATRIX 

49 DO 53 1=1, N 

53 WRI TE ( 5 , 103 ) (Q(I,J),J=1,N) 

C COMPUTE THE MATRIX QTRANSP0SE*0 

CALL MAT A ( Q ,GC »N ,N, 0 ) ■ 

C GET THE GENERAL STORAGE MODE FOR QTRANSPOSE*Q 

CALL MS TR ( GC . 0 « N* 1 . 6 ) 

WRITE(5,101) 

WRITE(5,212) 

212 FORMAT ( ' U TRANSPOSE Q IS’,/) 

C WRITE QTRANSPOSE+Q 

DO 51 1=1, N 

51 WR I TE ( 5 , 103 ) ( Q ( I , J ) , J= 1 , N ) 

WR I TE ( 5 , 10 1 ) 

C COMPUTE THE EIGENVALUES OF QTRANSPOSE*Q 

CALL E I GEN ( GC »D»N, 1 ) 

C COMPUTE THE ESTIMATE OF THE NORM OF THE MATRIX U 

E ( K ) =SQRT ( GC ( 1 ) ) 

WRIT E ( 5,214) E(K) 

214 FORMAT (• THE ESTIMATE OF THE NORM OF THE MATRIX IS'»F8.4) 

52 CONTINUE 

AGGREGATION 
WRITE(5,215) 

215 FORMAT ( ' 1 THE AGGREGATION MATRIX A AS A FUNCTION OF COUPLING PARAME 
1TAR ZETA IS',//) 

A 1 1 = -E TA3 ( 1 ) /ETA2< 1 ) 

A 1 =E ( l )*ETA4(1 ) /ETA1 ( 1 ) 

A12=E(2)*ETA4(l)/ETAl(2) 

A21 = E(4)*.ETA4(2)/ETA1 ( 1) 

A22 = -ETA3(2 )/ETA2(2 ) 

A2=E(3)*ETA4(2)/ETA1(2) 
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PASSIVE STABILITY OF THE SPINNING SKYLAB 
WRITE (5, 106) All ,A1,A12 

106 FORMAT ( F 12 • 6 * ' + • , F8 . 2 , • *Z ET A*Z ET A • , 20X , F8 . 2 , • * Z ETA • ) 

WR I TE ( 5 » 1 0 1 ) 

WRITE(5»107) A21,A22»A2 

107 FORMAT ( F8 . 2 » ' *ZET A ' ,20X,F12.6» • ♦ • , F 8 . 2 , • *Z ET A*Z ET A • ) 

ZETA* 1 . OE-7 

A(1,1)=A11+ZETA*ZETA*A1 

A( 1,2)=A12*ZETA 

A(2, 1 )=A21*ZETA 

A(2,2)=A22+ZETA*ZETA*A2 

DETA = A( 1 ,1 )*A(2,2)-A( 1 ,2 )*A(2, 1) 

WRITE (5,102) 

WRITE (5,220) ZETA 

220 FORMAT ( * THE AGGREGATION MATRIX A, FOR ZETA*' , EiO.4, ' IS* ,//) 

DO 60 1=1,2 

60 WRITE(5,221) ( A ( I » J ) » J= 1 » 2 ) 

221 FORMAT ( E 1 2 • 4 , 10X, E 1 2 • 4 , / ) 

WRITE(5,102) 

WRITE(5,222) A(1,1),DETA. 

222 FORMAT ( * A(l,l) IS ' , E 12 . 4 , ' , ANC DETERMINANT OF A IS ',E10.4,//) 
I F ( A ( 1 , 1 ) ) 61,62,62 

61 IF(DETA) 62,62,64 
64 WRITE(5,223) 

223 FORMAT { * THE OVERALL SYSTEM IS ASYMPTOT I CALLY STABLE ' ) 

GO TO 15 

62 WR I TE ( 5 , 224 ) 

224 FORMAT ( • LOOK FOR SMALLER VALUE OF THE COUPLING PARAMETAR ZETA') 

15 CALL EXIT : ‘ ' 1 * >! ‘ ' '' • 

END 

FEATURES SUPPORTED 

ONE WORD INTEGERS . 

IOCS 

CORE REQUIREMENTS FOR 

COMMON 0 VARIABLES 410 PROGRAM 2238 

END. OF COMPILATION . - . 

// XEQ.,- ; ....... 
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ESTIMATING THE NORMS OF THE INTERCONNECTING MATRICES 
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THE AGGREGATION MATRIX A AS A FUNCTION OF COUPLING PARAMETAR 2ETA IS 
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